The scripts can be found in the following repository

Link to the function: https://github.com/Metacommunity-Lab/HYD_UK-RI_FreswHabitats

Link to supplementary material: https://doi.org/10.1007/s10750-022-05123-0

Note to the version: All the .RData files are too heavy to be uploaded into the GitHub repository and can be downloaded from the following zenodo repository: https://zenodo.org/record/8155013

# Load basic functions to run all the simulations
library(dplyr); library(tibble); library(ggplot2); library(gridExtra);library(geosphere);library(sp);library(stars)
load("Paisaje_UK.RData")

The following tutorial is structured in 3 main parts corresponding to the main steps of the analysis:

1. Charging and preparation of the datasets. The charged datasets correspond to the different water covers (Permannet, Temporal, Ephemeral) and the corresponding databases matching among them. Distances calculations and network (graph) creation. All these initial steps were meant to prepare the database and explore the distribution of the habitats through the United kingdom and the Republic of Ireland.

2. Running the coalescent models.
   We run these models for each one of the explored combinations of habitats:

• Permanent without considering Temporal and Ephemeral
• Permanent considering Temporal and Ephemeral
• Temporal without considering Permanent and Ephemeral
• Temporal considering Permanent and Ephemeral
• Ephemeral without considering Temporal and Permanent
• Ephemeral considering Temporal and Permanent For each one of the runs of the current models we obtained the expected mean richness that later ploted and compared along figures.

3. Plotting and comparisons of the results. Plots scripts contain also the analysis carried for each comparison (Figure 5) where we crossed the richness values obtained from considering or not the other types of environments when assessing richness values for a determined habitat.

All freshwaters

We build and charge all the datasets needed to run the current model and run the first iteration model considering all the types of habitats together as well as the corresponding Filters previously calculated.

M.dist.UK<-as.matrix(dist(M.UK[,2:3], method = "euclidean"))

Spp.pool<-rep(1,200) # Uniform species abundance distribution is assumed for the pool
mod.temp<-rep(1,nrow(M.UK))
mod.temp[1:10]<-2
# First iteration model considering all freshwaters
## Filter comes from H2020_UK&IR_Fitness_set_up.R
Uk.metacomm<-iteration.model(replicas=14, Meta.pool=Spp.pool, m.pool=0.01, Js=M.UK[,6], id.module=M.UK[,7], filter.env=Filter,
                          M.dist=M.dist.UK, D50=4, m.max=1,
                          id.fixed=NULL, D50.fixed=0, m.max.fixed=0, comm.fixed=Spp.pool,
                          Lottery=T, it=100, prop.dead.by.it=0.05, id.obs=1:nrow(M.UK))

UK.out<-resume.out(Uk.metacomm)
names(UK.out)
UK.out.median<-UK.out[[1]]
N<-nrow(M.dist.UK)
S<-UK.out.median[11:(10+N)]
tempo<-cbind(M.UK.cent.D50.4[id.0,],0)
colnames(tempo)[ncol(tempo)]<-"S"
M.UK.S<-cbind(M.UK,S)
M.UK.S<-rbind(M.UK.S,tempo)
which(M.UK.S$grado.in.D50.4==min(M.UK.S$grado.in.D50.4))
M.UK.S<-M.UK.S[-2191,]
which(M.UK.S$grado.in.D50.4==min(M.UK.S$grado.in.D50.4))
M.UK.S<-M.UK.S[-5585,]
which(M.UK.S$grado.in.D50.4==min(M.UK.S$grado.in.D50.4))
M.UK.S<-M.UK.S[-4462,]

id.1<-which(M.UK[,7]==1)
id.2<-which(M.UK[,7]==2)
id.3<-which(M.UK[,7]==3)

id.1<-which(M.UK[,7]==1)
id.2<-which(M.UK[,7]==2)
id.3<-which(M.UK[,7]==3)
id_without_Ephemeral<-c(id.2,id.3)
id_without_temporal<-c(id.1,id.3)
id_without_permanent<-c(id.1,id.2)
id_permanet<-c(id.3)
id_temporal<-c(id.2)
id_ephemeral<-c(id.1)

Ephemeral freshwaters

We run again the current model to obtain the richness patterns only considering the Ephemeral freshwaters.

# Second iteration model considering only ephemeral freshwaters
## Filter comes from H2020_UK&IR_Fitness_set_up.R
Uk.metacomm_without_only_ephemeral<-iteration.model(replicas=14, Meta.pool=Spp.pool, m.pool=0.01, Js=M.UK[id_ephemeral,6], 
                                               id.module=M.UK[id_ephemeral,7], filter.env=Filter[,id_ephemeral],
                                               M.dist=M.dist.UK[id_ephemeral,id_ephemeral], D50=4, m.max=1,
                                               id.fixed=NULL, D50.fixed=0, m.max.fixed=0, comm.fixed=Spp.pool,
                                               Lottery=T, it=100, prop.dead.by.it=0.05, id.obs=1:length(id_ephemeral))
Uk.metacomm_without_only_ephemeral.out<-resume.out(Uk.metacomm_without_only_ephemeral)
S_ephy<-Uk.metacomm_without_only_ephemeral.out[[1]][10:(9+length(id_ephemeral))]
M.UK.S<-cbind(M.UK.S,0)
M.UK.S[id_ephemeral,ncol(M.UK.S)]<-S_ephy
colnames(M.UK.S)[ncol(M.UK.S)]<-"S_only_ephemeral"

Temporal freshwaters

We run again the current model to obtain the richness patterns only considering the Temporal freshwaters.

# Third iteration model considering only ephemeral freshwaters
## Filter comes from H2020_UK&IR_Fitness_set_up.R
Uk.metacomm_without_only_temporal<-iteration.model(replicas=14, Meta.pool=Spp.pool, m.pool=0.01, Js=M.UK[id_temporal,6], 
                                                    id.module=M.UK[id_temporal,7], filter.env=Filter[,id_temporal],
                                                    M.dist=M.dist.UK[id_temporal,id_temporal], D50=4, m.max=1,
                                                    id.fixed=NULL, D50.fixed=0, m.max.fixed=0, comm.fixed=Spp.pool,
                                                    Lottery=T, it=100, prop.dead.by.it=0.05, id.obs=1:length(id_temporal))
Uk.metacomm_without_only_temporal.out<-resume.out(Uk.metacomm_without_only_temporal)
S_temporal<-Uk.metacomm_without_only_temporal.out[[1]][10:(9+length(id_temporal))]
M.UK.S<-cbind(M.UK.S,0)
M.UK.S[id_temporal,ncol(M.UK.S)]<-S_temporal
colnames(M.UK.S)[ncol(M.UK.S)]<-"S_only_temporal"

Permanent freshwaters

We run again the current model to obtain the richness patterns only considering the Permanent freshwaters.

# Fourth iteration model considering only ephemeral freshwaters
## Filter comes from H2020_UK&IR_Fitness_set_up.R
Uk.metacomm_without_only_permanent<-iteration.model(replicas=14, Meta.pool=Spp.pool, m.pool=0.01, Js=M.UK[id_permanet,6], 
                                                   id.module=M.UK[id_permanet,7], filter.env=Filter[,id_permanet],
                                                   M.dist=M.dist.UK[id_permanet,id_permanet], D50=4, m.max=1,
                                                   id.fixed=NULL, D50.fixed=0, m.max.fixed=0, comm.fixed=Spp.pool,
                                                   Lottery=T, it=100, prop.dead.by.it=0.05, id.obs=1:length(id_permanet))

Uk.metacomm_without_only_permanent.out<-resume.out(Uk.metacomm_without_only_permanent)
S_permanent<-Uk.metacomm_without_only_permanent.out[[1]][10:(9+length(id_permanet))]
M.UK.S<-cbind(M.UK.S,0)
M.UK.S[id_permanet,ncol(M.UK.S)]<-S_permanent
colnames(M.UK.S)[ncol(M.UK.S)]<-"S_only_permanent"

Figure 1

Simulated species performances along ephemeral, temporal, and permanent freshwaters. The x-axis is the species ID—an identification number—and the y-axis is the species performance. Performance bias indicated the recruitment of individuals of each species in each freshwater environment in the coalescent assembly of communities.

# Figure 1 Fitness curves ####
S<-200       # Number of species
b1=0.05; a1<- -S*0.65*b1
Species<-1:S
Perf.1<- exp(a1+b1*Species)/(1+exp(a1+b1*Species)) # Performance of species pool 1

b2= -0.05; a2<- -S*0.35*b2
Perf.2<- exp(a2+b2*Species)/(1+exp(a2+b2*Species))# performance of species pool 2

Perf.matrix<-cbind(Perf.1, Perf.2)
colnames(Perf.matrix)<-c(20,1)

y<-c(rep(0,19),1,1,rep(0,19),rep(1,58),rep(0,4),rep(1,58),rep(0,19),1,1,rep(0,19))
y[9]<-1;y[187]<-1 #;y[99:101]<-0;
summary(glm(y~Species+I(Species^2), family = binomial))
p<-coefficients(glm(y~Species+I(Species^2), family = binomial))
Perf.3<- exp(p[1]+p[2]*Species+p[3]*Species^2)/(1+exp(p[1]+p[2]*Species+p[3]*Species^2))# performance species pool 3

#Plot
png(filename ="Figure1_Performance.png", 
    width = 700*3, height = 500*3, 
    units = "px",res = 300)
data.frame(Species,Perf.1,Perf.2,Perf.3)%>%
  ggplot()+
  geom_line(aes(x=Species, y=Perf.1), size=3,colour=cividis(3)[1],alpha=1)+
  geom_line(aes(x=Species, y=Perf.2), size=3,colour=cividis(3)[2],alpha=1)+
  geom_line(aes(x=Species, y=Perf.3), size=3,colour=cividis(3)[3],alpha=1)+
  scale_y_continuous(expand = c(0.02,0.02),breaks = c(0,0.25,0.5,0.75,1),labels = c(0,0.25,0.5,0.75,1))+
  scale_x_continuous(expand = c(0.01,0.01),breaks = c(1,25,50,75,100,125,150,175,200),labels = c(1,25,50,75,100,125,150,175,200))+
  labs(title = "Species habitat performance", y="Affinity")+
  annotate(x=25,y=1, label="Ephemeral",geom="label", fill=cividis(3)[2], colour="white",alpha=1, size=5)+
  annotate(x=100,y=0.9,label="Temporal", geom="label", fill=cividis(3)[3], colour="white",alpha=1, size=5)+
  annotate(x=175,y=1,label="Permanent",geom="label", fill=cividis(3)[1], colour="white",alpha=1, size=5)+
  theme_classic()
dev.off()

Figure 2A

Expected richness in ephemeral (A), temporal (B), and permanent freshwaters (C) for the United Kingdom and Ireland in a grid cell (10 × 10 km). Richness was estimated with a coalescent model that considered dispersal only (i) through the same freshwater environment—S unique—or (ii) through the different kinds of freshwater environments (ephemeral, temporal and permanent)— S all . The log ratio between these two estimations (iii)—S LR—indicates the relative effect of dispersal between freshwater environments on diversity.

#Ephemeral waters with temporal and permanent water dispersal
# Small histogram located at the top ____
small_histo<-ggplot(data.frame(M.UK.S[id_ephemeral,]))+
             aes(S)+
             geom_histogram(col="grey70", fill="grey70", alpha=1,binwidth = 1)+
             labs(x="Local Richness", y="counts")+
             #geom_density(aes(y =..density..* (nrow(data.frame(M.UK.S[id_ephemeral,])))),alpha = 0.1, fill = "navy")+
             theme_classic()

# PLot Ephemeral with the others
A <- ggplot() + geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black')+ 
  coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59))+
  theme_void()+
  geom_point(data = data.frame(M.UK.S[id_ephemeral,]), 
             aes(x = CENTROID_X, y =CENTROID_Y, color = (S_ephy)),shape=15, size=1.25,alpha=0.6) +
  labs(title="A)", subtitle = "Ephemeral without the others")+
  scale_color_viridis(name="S. unique Eph.",limits=c(0,35))+
  #scale_color_CUNILLERA(palette = "wildfire",discrete = F,name="S. unique Eph.")+ 
  annotation_custom(ggplotGrob(small_histo),  xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)


#Ephemeral waters without the others ____
# Small histogram located at the top
small_histo<-ggplot(data.frame(M.UK.S[id_ephemeral,]))+
             aes(S_ephy)+
             geom_histogram(col="grey70", fill="grey70", alpha=1, binwidth = 1)+
             labs(x="Local Richness", y="counts")+
            #geom_density(aes(y =..density..* (nrow(data.frame(M.UK.S[id_ephemeral,])))),alpha = 0.1, fill = "navy")+
             theme_classic()

# PLot Ephemeral alone
B <-ggplot()+
     geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black') + 
     coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59)) +theme_void()+
     geom_point(data = data.frame(M.UK.S[id_ephemeral,]),aes(x = CENTROID_X, y =CENTROID_Y,
                                                             color = (S)),shape=15, size=1.25,alpha=0.6)+ 
     labs(title="B)", subtitle = "Ephemeral with Temporal and Permanent")+
     scale_color_viridis(name="S. all Eph.", limits=c(0,35))+
     #scale_color_CUNILLERA(palette = "wildfire",discrete = F, name="S. all Eph.")+  
     #scale_color_gradient(low="yellow",high="red", guide_colourbar(title ="S"))+
     annotation_custom(ggplotGrob(small_histo),  xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)


# Difference in richness ____
# Data calculation
M.UK.S[,13]<-0
M.UK.S[id_ephemeral,13] <- log10((M.UK.S$S[id_ephemeral]/S_ephy))
colnames(M.UK.S)[13]<-"Diff.Ephe.log.ratio"

# Small histogram located at the top
small_histo<-ggplot(data.frame(M.UK.S[id_ephemeral,]))+aes(Diff.Ephe.log.ratio)+
               geom_histogram(col="grey70", fill="grey70", alpha=1, bins = 15)+
               labs(x="Richness ratio", y="counts")+
               #geom_density(aes(y =..density..*120),alpha = 0.1, fill = "navy")+
               theme_classic()



# PLot Logratio of the difference S/Seph
C <-ggplot()+
    geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black') + 
    coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59)) +theme_void()+
    geom_point(data = data.frame(M.UK.S[id_ephemeral,]), 
    aes(x = CENTROID_X, y =CENTROID_Y, color = (Diff.Ephe.log.ratio)),shape=15, size=1.25,alpha=0.6) + 
    labs(title="C)", subtitle = "Richness log ratio")+
    scale_color_gradient2(midpoint = 0,
                          low=viridis(1,direction = -1),
                          mid="white",
                          high=viridis(1,direction = 1),
                          guide_colourbar(title ="SLR"))+
    annotation_custom(ggplotGrob(small_histo), xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)


png(filename ="Figure2_Ephemeral.png", 
    width = 900*4.5, height = 600*4.5, 
    units = "px",res = 300)
grid.arrange(A,B,C, ncol=3, nrow=1)
dev.off()

Figure 2B

Expected richness in ephemeral (A), temporal (B), and permanent freshwaters (C) for the United Kingdom and Ireland in a grid cell (10 × 10 km). Richness was estimated with a coalescent model that considered dispersal only (i) through the same freshwater environment—S unique—or (ii) through the different kinds of freshwater environments (ephemeral, temporal and permanent)— S all . The log ratio between these two estimations (iii)—S LR—indicates the relative effect of dispersal between freshwater environments on diversity.

#Temporal waters with temporal and permanent water dispersal
# Small histogram located at the top ____
small_histo<-ggplot(data.frame(M.UK.S[id_temporal,]))+
  aes(S)+
  geom_histogram(col="grey70", fill="grey70", alpha=1,binwidth = 1)+
  labs(x="Local Richness", y="counts")+
  #geom_density(aes(y =..density..* (nrow(data.frame(M.UK.S[id_temporal,])))),alpha = 0.1, fill = "navy")+
  theme_classic()

# PLot Ephemeral with the others
A <- ggplot() + geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black')+ 
  coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59))+
  theme_void()+
  geom_point(data = data.frame(M.UK.S[id_temporal,]), 
             aes(x = CENTROID_X, y =CENTROID_Y, color = (S_temporal)),shape=15, size=1.25,alpha=0.6) +
  labs(title="A)", subtitle = "Temporal without the others")+
  scale_color_viridis(name="S. unique Temp.",limits=c(0,55))+
  #scale_color_CUNILLERA(palette = "wildfire",discrete = F,name="S. unique Temp.")+ 
  annotation_custom(ggplotGrob(small_histo),  xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)

  #Ephemeral waters without the others ____
  # Small histogram located at the top
small_histo<-ggplot(data.frame(M.UK.S[id_temporal,]))+
  aes(S_temporal)+
  geom_histogram(col="grey70", fill="grey70", alpha=1, binwidth = 1)+
  labs(x="Local Richness", y="counts")+
  #geom_density(aes(y =..density..* (nrow(data.frame(M.UK.S[id_temporal,])))),alpha = 0.1, fill = "navy")+
  theme_classic()

# PLot Ephemeral alone
B <-ggplot()+
  geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black') + 
  coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59)) +theme_void()+
  geom_point(data = data.frame(M.UK.S[id_temporal,]),aes(x = CENTROID_X, y =CENTROID_Y,
                                                          color = (S)),shape=15, size=1.25,alpha=0.6)+ 
  labs(title="B)", subtitle = "Temporal with Ephemeral and Permanent")+
  scale_color_viridis(name="S. all Temp.",limits=c(0,55))+
  #scale_color_CUNILLERA(palette = "wildfire",discrete = F, name="S. all Temp.")+  
  #scale_color_gradient(low="yellow",high="red", guide_colourbar(title ="S"))+
  annotation_custom(ggplotGrob(small_histo),  xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)

# Difference in richness ____
# Data calculation
M.UK.S[,13]<-0
M.UK.S[id_temporal,13] <- log10((M.UK.S$S[id_temporal]/S_temporal))
colnames(M.UK.S)[13]<-"Diff.Temp.log.ratio"

# Small histogram located at the top
small_histo<-ggplot(data.frame(M.UK.S[id_temporal,]))+aes(Diff.Temp.log.ratio)+
  geom_histogram(col="grey70", fill="grey70", alpha=1, bins = 15)+
  labs(x="Richness ratio", y="counts")+
  #geom_density(aes(y =..density..*200),alpha = 0.1, fill = "navy")+
  theme_classic()

# PLot Logratio of the difference S/Seph
C <-ggplot()+
  geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black') + 
  coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59)) +theme_void()+
  geom_point(data = data.frame(M.UK.S[id_temporal,]), 
             aes(x = CENTROID_X, y =CENTROID_Y, color = (Diff.Temp.log.ratio)),shape=15, size=1.25,alpha=0.6) + 
  labs(title="C)", subtitle = "Richness log ratio")+
  scale_color_gradient2(midpoint = 0,
                        low=viridis(1,direction = -1),
                        mid="white",
                        high=viridis(1,direction = 1),
                        guide_colourbar(title ="SLR"))+
  annotation_custom(ggplotGrob(small_histo), xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)



png(filename ="Figure3_Temporal.png", 
    width = 900*4.5, height = 600*4.5, 
    units = "px",res = 300)
grid.arrange(A,B,C, ncol=3, nrow=1)
dev.off()

Figure 2C

Expected richness in ephemeral (A), temporal (B), and permanent freshwaters (C) for the United Kingdom and Ireland in a grid cell (10 × 10 km). Richness was estimated with a coalescent model that considered dispersal only (i) through the same freshwater environment—S unique—or (ii) through the different kinds of freshwater environments (ephemeral, temporal and permanent)— S all . The log ratio between these two estimations (iii)—S LR—indicates the relative effect of dispersal between freshwater environments on diversity.

#Ephemeral waters with temporal and permanent water dispersal
# Small histogram located at the top ____
small_histo<-ggplot(data.frame(M.UK.S[id_permanet,]))+
  aes(S)+
  geom_histogram(col="grey70", fill="grey70", alpha=1,binwidth = 1)+
  labs(x="Local Richness", y="counts")+
  #geom_density(aes(y =..density..* (nrow(data.frame(M.UK.S[id_permanet,])))),alpha = 0.1, fill = "navy")+
  theme_classic()

# PLot Ephemeral with the others
A <-ggplot() + geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black')+ 
  coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59))+
  theme_void()+
  geom_point(data = data.frame(M.UK.S[id_permanet,]), 
             aes(x = CENTROID_X, y =CENTROID_Y, color = (S_permanent)),shape=15, size=1.25,alpha=0.6) +
  labs(title="A)", subtitle = "Permanent without the others")+
  scale_color_viridis(name="S. unique Perm.",limits=c(0,80))+
  #scale_color_CUNILLERA(palette = "wildfire",discrete = F,name="S. unique Perm.")+ 
  annotation_custom(ggplotGrob(small_histo),  xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2) 

  #Ephemeral waters without the others ____
  # Small histogram located at the top
  small_histo<-ggplot(data.frame(M.UK.S[id_permanet,]))+
    aes(S_permanent)+
    geom_histogram(col="grey70", fill="grey70", alpha=1,binwidth = 1)+
    labs(x="Local Richness", y="counts")+
    #geom_density(aes(y =..density..* (nrow(data.frame(M.UK.S[id_permanet,])))),alpha = 0.1, fill = "navy")+
    theme_classic()


# PLot Ephemeral alone
B <-ggplot()+
  geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black') + 
  coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59)) +theme_void()+
  geom_point(data = data.frame(M.UK.S[id_permanet,]),aes(x = CENTROID_X, y =CENTROID_Y,
                                                          color = (S)),shape=15, size=1.25,alpha=0.6)+ 
  labs(title="B)", subtitle = "Permanent with Ephemeral and Temporal")+
  scale_color_viridis(name="S. all Perm.",limits=c(0,80))+
  #scale_color_CUNILLERA(palette = "wildfire",discrete = F, name="S. all Perm.")+  
  #scale_color_gradient(low="yellow",high="red", guide_colourbar(title ="S"))+
  annotation_custom(ggplotGrob(small_histo),  xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)


# Difference in richness ____
# Data calculation
M.UK.S[,13]<-0
M.UK.S[id_permanet,13] <- log10((M.UK.S$S[id_permanet]/S_permanent))
colnames(M.UK.S)[13]<-"Diff.Perm.log.ratio"

# Small histogram located at the top
small_histo<-ggplot(data.frame(M.UK.S[id_permanet,]))+aes(Diff.Perm.log.ratio)+
  geom_histogram(col="grey70", fill="grey70", alpha=1,bins = 15)+
  labs(x="Richness ratio", y="counts")+
  #geom_density(aes(y =..density..*90),alpha = 0.1, fill = "navy")+
  theme_classic()

# PLot Logratio of the difference S/Seph
C <-ggplot()+
  geom_polygon(data = worldmap, aes(x = long, y = lat, group = group),fill = 'gray90', color = 'black') + 
  coord_fixed(ratio = 1.5, xlim = c(-10.5,2), ylim = c(50, 59)) +theme_void()+
  geom_point(data = data.frame(M.UK.S[id_permanet,]), 
             aes(x = CENTROID_X, y =CENTROID_Y, color = (Diff.Perm.log.ratio)),shape=15, size=1.25,alpha=0.6) + 
  labs(title="C)", subtitle = "Richness log ratio")+
  scale_color_gradient2(midpoint = 0,
                        low=viridis(1,direction = -1),
                        mid="white",
                        high=viridis(1,direction = 1),
                        guide_colourbar(title ="SLR"))+
  annotation_custom(ggplotGrob(small_histo), xmin = -1.5, xmax = 4.2,  ymin = 56.1, ymax = 59.2)



png(filename ="Figure4_Permanent.png", 
    width = 900*4.5, height = 600*4.5, 
    units = "px",res = 300)
grid.arrange(A,B,C, ncol=3, nrow=1)
dev.off()

Figure 3

Effect of the connection among freshwater environments on local diversity as a function of local diversity when this connection does not occur. The red dashed line indicates the 1:1 relationship, e.g., no effect of interfreshwater dispersal. Scatterplots show that the coupling between different types of environments foster local diversity in a magnitude that increases from ephemeral to permanent waters and from poorer to richer communities.

# Scatterplot S/Sephy ____
D_eph <- ggplot(data.frame(M.UK.S[id_ephemeral,]))+ aes(x=S_ephy, y=S)+
  geom_point(aes(fill=S),shape=22,size=3, col="grey20", alpha=0.2)+
  geom_smooth(method="lm", color="black", size=1.5)+
  geom_abline(intercept=0, slope=1, color="red", size=1, linetype=2)+
  scale_fill_viridis()+
  #scale_fill_CUNILLERA(palette = "wildfire",discrete = F)+
  labs(title="A)", subtitle = "Richness relationship (S. all Eph./S. unique Eph.)",
       y="S. all Eph.",x="S. unique Eph.")+
  geom_text(aes(y=40, x=10), label=paste("Slope=", 
                                        as.character(round(lm(M.UK.S[id_ephemeral,]$S~S_ephy)$coefficients[2],2)),
                                        sep=" "))+
  scale_y_continuous(limits = c(0,40))+
  scale_x_continuous(limits = c(0,40))+
  theme_classic()+
  theme(plot.margin=margin(t = 2, r = 2.5, b = 0.5, l = 2.5, unit = "cm"),
        legend.position = "none")

# Scatterplot S/Sephy ____
D_temp <- ggplot(data.frame(M.UK.S[id_temporal,]))+ aes(x=S_temporal, y=S)+
  geom_point(aes(fill=S),shape=22,size=3, col="grey20", alpha=0.2)+
  geom_smooth(method="lm", color="black", size=1.5)+
  geom_abline(intercept=0, slope=1, color="red", size=1, linetype=2)+
  scale_fill_viridis()+
  #scale_fill_CUNILLERA(palette = "wildfire",discrete = F)+
  labs(title="B)", subtitle = "Richness relationship (S. all Temp./S. unique Temp.)",
       y="S. all Temp.",x="S. unique Temp.")+
  geom_text(aes(y=60, x=10), label=paste("Slope=", 
                                        as.character(round(lm(M.UK.S[id_temporal,]$S~S_temporal)$coefficients[2],2)),
                                        sep=" "))+
  scale_y_continuous(limits = c(0,60))+
  scale_x_continuous(limits = c(0,60))+
  theme_classic()+
  theme(plot.margin=margin(t = 2, r = 2.5, b = 0.5, l = 2.5, unit = "cm"),
        legend.position = "none")

# Scatterplot S/Sperm ____
D_perm <- ggplot(data.frame(M.UK.S[id_permanet,]))+ aes(x=S_permanent, y=S)+
  geom_point(aes(fill=S),shape=22,size=3, col="grey20", alpha=0.2)+
  geom_smooth(method="lm", color="black", size=1.5)+
  geom_abline(intercept=0, slope=1, color="red", size=1, linetype=2)+
  scale_fill_viridis()+
  #scale_fill_CUNILLERA(palette = "wildfire",discrete = F)+
  labs(title="C)", subtitle = "Richness relationship (S. all Perm./S. unique Perm.)",
       y="S. all Perm.",x="S. unique Perm.")+
  geom_text(aes(y=100, x=10), label=paste("Slope=", 
                                        as.character(round(lm(M.UK.S[id_permanet,]$S~S_permanent)$coefficients[2],2)),
                                        sep=" "))+
  scale_y_continuous(limits = c(0,100))+
  scale_x_continuous(limits = c(0,100))+
  theme_classic()+
  theme(plot.margin=margin(t = 2, r = 2.5, b = 0.5, l = 2.5, unit = "cm"),
        legend.position = "none")

png(filename ="Figure5_Rel.png", 
    width = 1000*7, height = 300*4.5, 
    units = "px",res = 300)
grid.arrange(D_eph,D_temp,D_perm, ncol=3, nrow=1)
dev.off()

table_result <- rbind(
as.matrix(round(summary(lm(M.UK.S[id_ephemeral,]$S~S_ephy))[[4]],3)),
as.matrix(round(summary(lm(M.UK.S[id_temporal,]$S~S_temporal))[[4]],3)),
as.matrix(round(summary(lm(M.UK.S[id_permanet,]$S~S_permanent))[[4]],3))
)

Supplementary Figure 1

Scenarios richness relationships along latitude (left column) and longitude (right column) for each aquatic habitat type: Ephemeral, Temporal, Permanent. Colors differentiate the two scenarios: S. unique where the model is run over a landscape with only one type of habitat and S. all where the model is run considering the three habitats (i.e., whole waterscape).

# Comparison between lat and long 

Plot_geogr_relationships <- gridExtra::grid.arrange(
# Ephemeral
gridExtra::arrangeGrob(
data.frame(lat=data.frame(M.UK.S[id_ephemeral,])[,5],long=data.frame(M.UK.S[id_ephemeral,])[,4],
           S=data.frame(M.UK.S[id_ephemeral,])[,12],Type="ephem_All") %>% 
bind_rows(data.frame(lat=data.frame(M.UK.S[id_ephemeral,])[,5],long=data.frame(M.UK.S[id_ephemeral,])[,4],
          S=S_ephy,Type="ephem_Only"))%>% 
  ggplot()+
  geom_jitter(aes(x=lat, y=S, fill=Type), shape=21, alpha=0.1, colour="black")+
  geom_density2d(aes(x=lat, y=S, colour=Type),alpha=.5, size=1.5, linetype=2)+ 
  geom_smooth(aes(x=lat, y=S, colour=Type), method = "lm")+
  scale_fill_viridis(discrete = T)+
  scale_colour_viridis(discrete = T, label=c("S. all Eph.", "S. unique Eph."))+
  guides(fill="none")+labs(x="Latitude", y="Species richness")+
  geom_text(aes(x=59,y=33), colour=viridis(1,direction = -1),
            label=paste("Slope=",
            round(as.numeric(lm(S_ephy~
            data.frame(M.UK.S[id_ephemeral,])[,5])[[1]][2]),2),
            sep=" "))+
  geom_text(aes(x=59,y=35), colour=viridis(1,direction = 1),
            label=paste("Slope=",
            round(as.numeric(lm(data.frame(M.UK.S[id_ephemeral,])[,12]~
            data.frame(M.UK.S[id_ephemeral,])[,5])[[1]][2]),2),
            sep=" "))+
  theme_classic(),

data.frame(lat=data.frame(M.UK.S[id_ephemeral,])[,5],long=data.frame(M.UK.S[id_ephemeral,])[,4],
           S=data.frame(M.UK.S[id_ephemeral,])[,12],Type="ephem_All") %>% 
  bind_rows(data.frame(lat=data.frame(M.UK.S[id_ephemeral,])[,5],long=data.frame(M.UK.S[id_ephemeral,])[,4],
                       S=S_ephy,Type="ephem_Only"))%>% 
  ggplot()+
  geom_jitter(aes(x=long, y=S, fill=Type), shape=21, alpha=0.1, colour="black")+
  geom_density2d(aes(x=long, y=S, colour=Type),alpha=.3, size=1.5, linetype=2)+ 
  geom_smooth(aes(x=long, y=S, colour=Type), method = "lm", size=2, linetype=1)+
  scale_fill_viridis(discrete = T)+
  scale_colour_viridis(discrete = T, label=c("S. all Eph.", "S. unique Eph."))+
  guides(fill="none")+labs(x="Longitude", y="Species richness")+
  geom_text(aes(x=-1,y=33), colour=viridis(1,direction = -1),
            label=paste("Slope=",
            round(as.numeric(lm(S_ephy~
            data.frame(M.UK.S[id_ephemeral,])[,4])[[1]][2]),2),
            sep=" "))+
  geom_text(aes(x=-1,y=35), colour=viridis(1,direction = 1),
            label=paste("Slope=",
            round(as.numeric(lm(data.frame(M.UK.S[id_ephemeral,])[,12]~
            data.frame(M.UK.S[id_ephemeral,])[,4])[[1]][2]),2),
            sep=" "))+
  theme_classic(),
ncol=2, top="Ephemeral"),

# Temporal 
gridExtra::arrangeGrob(
  data.frame(lat=data.frame(M.UK.S[id_temporal,])[,5],long=data.frame(M.UK.S[id_temporal,])[,4],
             S=data.frame(M.UK.S[id_temporal,])[,12],Type="tempo_All") %>% 
    bind_rows(data.frame(lat=data.frame(M.UK.S[id_temporal,])[,5],long=data.frame(M.UK.S[id_temporal,])[,4],
                         S=S_temporal,Type="tempo_Only"))%>% 
    ggplot()+
    geom_jitter(aes(x=lat, y=S, fill=Type), shape=21, alpha=0.1, colour="black")+
    geom_density2d(aes(x=lat, y=S, colour=Type),alpha=.5, size=1.5, linetype=2)+ 
    geom_smooth(aes(x=lat, y=S, colour=Type), method = "lm")+
    scale_fill_viridis(discrete = T)+
    scale_colour_viridis(discrete = T, label=c("S. all Temp.", "S. unique Temp."))+
    guides(fill="none")+labs(x="Latitude", y="Species richness")+
    geom_text(aes(x=59,y=50), colour=viridis(1,direction = -1),
              label=paste("Slope=",
              round(as.numeric(lm(S_temporal~
              data.frame(M.UK.S[id_temporal,])[,5])[[1]][2]),2),
              sep=" "))+
    geom_text(aes(x=59,y=53), colour=viridis(1,direction = 1),
              label=paste("Slope=",
              round(as.numeric(lm(data.frame(M.UK.S[id_temporal,])[,12]~
              data.frame(M.UK.S[id_temporal,])[,5])[[1]][2]),2),
              sep=" "))+
    theme_classic(),
  
  data.frame(lat=data.frame(M.UK.S[id_temporal,])[,5],long=data.frame(M.UK.S[id_temporal,])[,4],
             S=data.frame(M.UK.S[id_temporal,])[,12],Type="tempo_All") %>% 
    bind_rows(data.frame(lat=data.frame(M.UK.S[id_temporal,])[,5],long=data.frame(M.UK.S[id_temporal,])[,4],
                         S=S_temporal,Type="tempo_Only"))%>% 
    ggplot()+
    geom_jitter(aes(x=long, y=S, fill=Type), shape=21, alpha=0.1, colour="black")+
    geom_density2d(aes(x=long, y=S, colour=Type),alpha=.3, size=1.5, linetype=2)+ 
    geom_smooth(aes(x=long, y=S, colour=Type), method = "lm", size=2, linetype=1)+
    scale_fill_viridis(discrete = T)+
    scale_colour_viridis(discrete = T, label=c("S. all Temp.", "S. unique Temp."))+
    guides(fill="none")+labs(x="Longitude", y="Species richness")+
    geom_text(aes(x=0,y=50), colour=viridis(1,direction = -1),
              label=paste("Slope=",
              round(as.numeric(lm(S_temporal~
              data.frame(M.UK.S[id_temporal,])[,4])[[1]][2]),2),
              sep=" "))+
    geom_text(aes(x=0,y=53), colour=viridis(1,direction = 1),
              label=paste("Slope=",
              round(as.numeric(lm(data.frame(M.UK.S[id_temporal,])[,12]~
              data.frame(M.UK.S[id_temporal,])[,4])[[1]][2]),2),
              sep=" "))+
    theme_classic(),
  ncol=2, top="Temporal"),

# Permanent
gridExtra::arrangeGrob(
  data.frame(lat=data.frame(M.UK.S[id_permanet,])[,5],long=data.frame(M.UK.S[id_permanet,])[,4],
             S=data.frame(M.UK.S[id_permanet,])[,12],Type="perma_All") %>% 
    bind_rows(data.frame(lat=data.frame(M.UK.S[id_permanet,])[,5],long=data.frame(M.UK.S[id_permanet,])[,4],
                         S=S_permanent,Type="perma_Only"))%>% 
    ggplot()+
    geom_jitter(aes(x=lat, y=S, fill=Type), shape=21, alpha=0.1, colour="black")+
    geom_density2d(aes(x=lat, y=S, colour=Type),alpha=.5, size=1.5, linetype=2)+ 
    geom_smooth(aes(x=lat, y=S, colour=Type), method = "lm")+
    scale_fill_viridis(discrete = T)+
    scale_colour_viridis(discrete = T, label=c("S. all Eph.", "S. unique Eph."))+
    guides(fill="none")+labs(x="Latitude", y="Species richness")+
    geom_text(aes(x=59,y=70), colour=viridis(1,direction = -1),
              label=paste("Slope=",
              round(as.numeric(lm(S_permanent~
              data.frame(M.UK.S[id_permanet,])[,5])[[1]][2]),2),
              sep=" "))+
    geom_text(aes(x=59,y=75), colour=viridis(1,direction = 1),
              label=paste("Slope=",
              round(as.numeric(lm(data.frame(M.UK.S[id_permanet,])[,12]~
              data.frame(M.UK.S[id_permanet,])[,5])[[1]][2]),2),
              sep=" "))+
    theme_classic(),
  
  data.frame(lat=data.frame(M.UK.S[id_permanet,])[,5],long=data.frame(M.UK.S[id_permanet,])[,4],
             S=data.frame(M.UK.S[id_permanet,])[,12],Type="perma_All") %>% 
    bind_rows(data.frame(lat=data.frame(M.UK.S[id_permanet,])[,5],long=data.frame(M.UK.S[id_permanet,])[,4],
                         S=S_permanent,Type="perma_Only"))%>% 
    ggplot()+
    geom_jitter(aes(x=long, y=S, fill=Type), shape=21, alpha=0.1, colour="black")+
    geom_density2d(aes(x=long, y=S, colour=Type),alpha=.3, size=1.5, linetype=2)+ 
    geom_smooth(aes(x=long, y=S, colour=Type), method = "lm", size=2, linetype=1)+
    scale_fill_viridis(discrete = T)+
    scale_colour_viridis(discrete = T, label=c("S. all Perm.", "S. unique Perm."))+
    guides(fill="none")+labs(x="Longitude", y="Species richness")+
    geom_text(aes(x=0,y=70), colour=viridis(1,direction = -1),
    label=paste("Slope=",
    round(as.numeric(lm(S_permanent~
                        data.frame(M.UK.S[id_permanet,])[,4])[[1]][2]),2),
    sep=" "))+
    geom_text(aes(x=0,y=75), colour=viridis(1,direction = 1),
    label=paste("Slope=",
    round(as.numeric(lm(data.frame(M.UK.S[id_permanet,])[,12]~
                        data.frame(M.UK.S[id_permanet,])[,4])[[1]][2]),2),
    sep=" "))+
    theme_classic(),
  ncol=2, top="Permanent"), 
nrow=3)

png(filename ="FigureSup_Rel.png", 
    width = 600*7, height = 600*7, 
    units = "px",res = 300)
grid.arrange(Plot_geogr_relationships)
dev.off()

Extra network plot

Network strcutre plots made from degree calculations to represent main centrality patterns accross the United Kingdom and Republic of Ireland. These plots were not included in the published ms but exemplify the distribution of the habitats and islands main centrality patterns for source, sink connectivity and stepping stone connectivity.

###########################################################################################
            ####### ####### ####### PLOTS ####### ####### #######
###########################################################################################
## M.final.UK esta es la matriz completa dierenciando cuerpos de agua y tierra
## M.final.FW.UK esta es la matriz de agua (toda junta) y tiera 
M.final.UK->MM
M.final.FW.UK->MM.FW


#### In_Degree_MassEffect #### D50=4Km
#mtext("In-degree (Mass effect) \n", outer = TRUE, cex = 1.2)
par(mar=c(2.5,2.5,1.5,1.5),bty="l",cex=1.8)
par(mfrow=c(3,4),oma = c(0, 0, 0, 0))
plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Ephemeral (4Km)")
ii <- cut(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==1)], 
          breaks = seq(min(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==1)]), 
                       max(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==1)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==1)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==1)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)
text(-8,60,"IN DEGREE\nMass effect", )
rm(ii)

plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Temporal (4Km)")
ii <- cut(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==2)], 
          breaks = seq(min(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==2)]), 
                       max(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==2)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==2)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==2)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)
rm(ii)

plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Permanent (4Km)")
ii <- cut(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==3)], 
          breaks = seq(min(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==3)]), 
                       max(MM$grado.in.D50.4[which(MM$`EF(1).T(2).P(3)`==3)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==3)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==3)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)

rm(ii)

plot(MM.FW$CENTROID_Y~MM.FW$CENTROID_X,col="gray",cex=0.4,pch=15,main="All-freshwater (4Km)")
ii <- cut(MM.FW$grado.in.D50.4, 
          breaks = seq(min(MM.FW$grado.in.D50.4), 
                       max(MM.FW$grado.in.D50.4), len = 600), include.lowest = TRUE)
points(MM.FW$CENTROID_Y~MM.FW$CENTROID_X,cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)


#### Out_Degree_SourcePatchEffect ####
rm(ii)
plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Ephemeral (4Km)")
ii <- cut(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==1)], 
          breaks = seq(min(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==1)]), 
                       max(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==1)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==1)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==1)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)
rm(ii)
text(-8,60,"OUT DEGREE\nsource patches \nin mass effect")


plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Temporal (4Km)")
ii <- cut(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==2)], 
          breaks = seq(min(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==2)]), 
                       max(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==2)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==2)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==2)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)
rm(ii)

plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Permanent (4Km)")
ii <- cut(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==3)], 
          breaks = seq(min(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==3)]), 
                       max(MM$grado.out.D50.4[which(MM$`EF(1).T(2).P(3)`==3)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==3)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==3)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)

rm(ii)

plot(MM.FW$CENTROID_Y~MM.FW$CENTROID_X,col="gray",cex=0.4,pch=15,main="All-freshwater (4Km)")
ii <- cut(MM.FW$grado.out.D50.4, 
          breaks = seq(min(MM.FW$grado.out.D50.4), 
                       max(MM.FW$grado.out.D50.4), len = 600), include.lowest = TRUE)
points(MM.FW$CENTROID_Y~MM.FW$CENTROID_X,cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)
#mtext("Out-degree (Source effect) \n", outer = TRUE, cex = 1.2)

##### bc
rm(ii)
plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Ephemeral (4Km)")
ii <- cut(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==1)], 
          breaks = seq(min(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==1)]), 
                       max(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==1)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==1)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==1)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)
rm(ii)
text(-10,60,"Betweenness\n steping stone \nlandscape dispersal", pos = 4)
plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Temporal (4Km)")
ii <- cut(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==2)], 
          breaks = seq(min(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==2)]), 
                       max(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==2)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==2)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==2)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)
rm(ii)

plot(MM$CENTROID_Y~MM$CENTROID_X,col="gray",cex=0.4,pch=15,main="Permanent (4Km)")
ii <- cut(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==3)], 
          breaks = seq(min(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==3)]), 
                       max(MM$bc.D50.4[which(MM$`EF(1).T(2).P(3)`==3)]), len = 600), include.lowest = TRUE)
points(MM$CENTROID_Y[which(MM$`EF(1).T(2).P(3)`==3)]~MM$CENTROID_X[which(MM$`EF(1).T(2).P(3)`==3)],cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)

rm(ii)

plot(MM.FW$CENTROID_Y~MM.FW$CENTROID_X,col="gray",cex=0.4,pch=15,main="All-freshwater (4Km)")
ii <- cut(MM.FW$bc.D50.4, 
          breaks = seq(min(MM.FW$bc.D50.4), 
                       max(MM.FW$bc.D50.4), len = 600), include.lowest = TRUE)
points(MM.FW$CENTROID_Y~MM.FW$CENTROID_X,cex=0.4,
       col = colorRampPalette(c("yellow", "orange", "red"))(599)[ii],pch=15)